<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  
  <title>twoPhaseEulerFoam 全解读之二 | Giskard&#39;s CFD Learning Tricks</title>
  <meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1">
  <meta name="description" content="本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇对 twoPhaseEulerFoam 中的 UEqn.H 和 pEqn.H 中的代码进行详细地的解读。">
<meta property="og:type" content="article">
<meta property="og:title" content="twoPhaseEulerFoam 全解读之二">
<meta property="og:url" content="http://xiaopingqiu.github.io/2015/05/17/twoPhaseEulerFoam2/index.html">
<meta property="og:site_name" content="Giskard's CFD Learning Tricks">
<meta property="og:description" content="本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇对 twoPhaseEulerFoam 中的 UEqn.H 和 pEqn.H 中的代码进行详细地的解读。">
<meta name="twitter:card" content="summary">
<meta name="twitter:title" content="twoPhaseEulerFoam 全解读之二">
<meta name="twitter:description" content="本系列将对OpenFOAM-2.1.1 中的 twoPhaseEulerFoam 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇对 twoPhaseEulerFoam 中的 UEqn.H 和 pEqn.H 中的代码进行详细地的解读。">
  
  
    <link rel="icon" href="/favicon.png">
  
  <link href="//fonts.googleapis.com/css?family=Source+Code+Pro" rel="stylesheet" type="text/css">
  <link rel="stylesheet" href="/css/style.css" type="text/css">
  
<!-- Google Analytics -->
<script type="text/javascript">
(function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
(i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
})(window,document,'script','//www.google-analytics.com/analytics.js','ga');

##ga('create', '[object Object]', 'auto');
ga('create', 'UA-62501686-1', 'auto');
ga('send', 'pageview');

</script>
<!-- End Google Analytics -->


</head>
<body>
  <div id="container">
    <div id="wrap">
      <header id="header">
  <div id="banner"></div>
  <div id="header-outer" class="outer">
    <div id="header-title" class="inner">
      <h1 id="logo-wrap">
        <a href="/" id="logo">Giskard&#39;s CFD Learning Tricks</a>
      </h1>
      
        <h2 id="subtitle-wrap">
          <a href="/" id="subtitle">CFD and Scientific Computing</a>
        </h2>
      
    </div>
    <div id="header-inner" class="inner">
      <nav id="main-nav">
        <a id="main-nav-toggle" class="nav-icon"></a>
        
          <a class="main-nav-link" href="/">Home</a>
        
          <a class="main-nav-link" href="/archives">Archives</a>
        
          <a class="main-nav-link" href="/atom.xml">Rss</a>
        
          <a class="main-nav-link" href="/about">About</a>
        
      </nav>
      <nav id="sub-nav">
        
        <a id="nav-search-btn" class="nav-icon" title="Search"></a>
      </nav>
      <div id="search-form-wrap">
        <form action="//google.com/search" method="get" accept-charset="UTF-8" class="search-form"><input type="search" name="q" results="0" class="search-form-input" placeholder="Search"><button type="submit" class="search-form-submit">&#xF002;</button><input type="hidden" name="q" value="site:http://xiaopingqiu.github.io"></form>
      </div>
    </div>
  </div>
</header>

      <div class="outer">
        <section id="main"><article id="post-twoPhaseEulerFoam2" class="article article-type-post" itemscope itemprop="blogPost">
  <div class="article-meta">
    <a href="/2015/05/17/twoPhaseEulerFoam2/" class="article-date">
  <time datetime="2015-05-17T07:23:04.000Z" itemprop="datePublished">2015-05-17</time>
</a>
    
  <div class="article-category">
    <a class="article-category-link" href="/categories/OpenFOAM/">OpenFOAM</a>
  </div>

  </div>
  <div class="article-inner">
    
    
      <header class="article-header">
        
  
    <h1 class="article-title" itemprop="name">
      twoPhaseEulerFoam 全解读之二
    </h1>
  

      </header>
    
    <div class="article-entry" itemprop="articleBody">
      
        <p>本系列将对OpenFOAM-2.1.1 中的 <code>twoPhaseEulerFoam</code> 求解器进行完全解读，共分三部分：方程推导，代码解读，补充说明。本篇对 <code>twoPhaseEulerFoam</code> 中的 <code>UEqn.H</code> 和 <code>pEqn.H</code> 中的代码进行详细地的解读。</p>
<a id="more"></a>
<h2 id="2-_代码解读">2. 代码解读</h2><h3 id="2-1-_UEqn">2.1. UEqn</h3><p>前一篇导出了分散相的动量守恒方程<br>$$<br>\begin{aligned}<br>&amp;(1+\frac{\alpha_b \rho_b}{\rho_a} C_{vm})(\frac{\partial U_a}{\partial t} + U_a\cdot \nabla U_a ) -\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ] + \nabla \cdot \left[ R_{c,a}\right] + \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ -\nu_{eff}\nabla U_a + R_{c,a}\right] \\<br>= &amp; -\frac{\alpha_b}{\rho_a} K U_a - \frac{\alpha_b}{\rho_a} \left\{ {C_l (\alpha_b \rho_b + \alpha_a \rho_a) U_r \times (\nabla \times U) -  C_{vm}\rho_b\left[ {\frac{\partial U_b}{\partial t} + U_b \cdot \nabla U_b }\right] } \right\} -\frac{\nabla p}{\rho_a} + g + \frac{\alpha_b}{\rho_a} K U_b<br> \end{aligned}<br>$$<br>这一篇分析<code>twoPhaseEulerFoam</code>求解器是怎么来对动量方程进行离散的，以及，如果通过构建压力方程来对速度进行修正以保证两相的连续性。</p>
<p>注意上述动量方程中，有两项还需要处理一下：<br>$$U_a \cdot \nabla U_a=\nabla\cdot(U_aU_a)-U_a(\nabla\cdot U_a)$$</p>
<p>$$\frac{\nabla(\alpha_a)}{\alpha_a}\cdot\left[-\nu_{eff} \nabla U_a\right] = \nabla \cdot\left[ -\nu_{eff}\frac{(\nabla \alpha_a)}{\alpha_a} U_a \right]-U_a\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right)$$<br>转化前后的形式从数学上来等价的，但是在有限体积离散过程中，转化前的 $U_a \cdot \nabla U_a$  和 $\frac{\nabla(\alpha_a)}{\alpha_a}\cdot\left[-\nu_{eff} \nabla U_a\right]$ 对于$U_a$来说是非守恒的，转化后的形式是守恒的(参考<a href="http://www.sciencedirect.com/science/article/pii/S0029549309003021" target="_blank" rel="external">这篇文献</a>，注意这样转化后，动量方程空间上是守恒的，但时间上仍是不守恒的，<a href="http://www.cfd-online.com/Forums/openfoam-solving/71141-rewriting-twophaseeulerfoam-conservative-form.html" target="_blank" rel="external">这个帖子</a>也有一些有价值的信息)。</p>
<p>这样转化以后，得到的动量方程就跟<code>twoPhaseEulerFoam</code>里的定义是一模一样了:<br>$$<br>\begin{aligned}<br>&amp;(1+\frac{\alpha_b \rho_b}{\rho_a} C_{vm})\left(\frac{\partial U_a}{\partial t} + \nabla\cdot(U_aU_a)-U_a(\nabla\cdot U_a) \right)\\<br>-&amp;\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ] + \nabla \cdot \left[ R_{c,a}\right] +\nabla \cdot\left[ -\nu_{eff}\frac{(\nabla \alpha_a)}{\alpha_a}(\nabla U_a)\right]-U_a\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right) \\<br>+ &amp; \frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ R_{c,a}\right] \\<br>= &amp; -\frac{\alpha_b}{\rho_a} K U_a - \frac{\alpha_b}{\rho_a} \left\{ {C_l (\alpha_b \rho_b + \alpha_a \rho_a) U_r \times (\nabla \times U) -  C_{vm}\rho_b\left[ {\frac{\partial U_b}{\partial t} + U_b \cdot \nabla U_b }\right] } \right\} \\<br>-&amp;\frac{\nabla p}{\rho_a} + g + \frac{\alpha_b}{\rho_a} K U_b<br> \end{aligned}<br>$$</p>
<p>下面将动量方程的每一项与<code>twoPhaseEulerFoam</code>的<code>UEqn.H</code>的代码一一对应。<br><figure class="highlight lisp"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br></pre></td><td class="code"><pre><span class="line"><span class="list">(<span class="keyword">scalar</span><span class="list">(<span class="number">1</span>)</span> + Cvm*rhob*beta/rhoa)</span>*</span><br><span class="line">            <span class="list">(</span><br><span class="line">                <span class="keyword">fvm</span>:<span class="keyword">:ddt</span><span class="list">(<span class="keyword">Ua</span>)</span></span><br><span class="line">              + fvm:<span class="keyword">:div</span><span class="list">(<span class="keyword">phia</span>, Ua, <span class="string">"div(phia,Ua)"</span>)</span></span><br><span class="line">              - fvm:<span class="keyword">:Sp</span><span class="list">(<span class="keyword">fvc</span>:<span class="keyword">:div</span><span class="list">(<span class="keyword">phia</span>)</span>, Ua)</span></span><br><span class="line">            )</span></span><br></pre></td></tr></table></figure></p>
<p><code>fvm::ddt(Ua)</code>对应$\frac{\partial U_a}{\partial t}$，<code>phia</code>定义为<code>fvc::interpolate(Ua) &amp; mesh.Sf()</code>，于是<code>fvm::div(phia, Ua, &quot;div(phia,Ua)&quot;)</code> 和 <code>fvm::Sp(fvc::div(phia), Ua)</code> 便分别对应 $\nabla\cdot(U_aU_a)$ 和 $U_a(\nabla\cdot U_a)$了 <strong>[ <em>注一</em> ]</strong>。</p>
<figure class="highlight stylus"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br></pre></td><td class="code"><pre><span class="line">- fvm::<span class="function"><span class="title">laplacian</span><span class="params">(nuEffa, Ua)</span></span></span><br><span class="line">+ fvc::<span class="function"><span class="title">div</span><span class="params">(Rca)</span></span></span><br><span class="line">+ fvm::<span class="function"><span class="title">div</span><span class="params">(phiRa, Ua, <span class="string">"div(phia,Ua)"</span>)</span></span></span><br><span class="line">- fvm::<span class="function"><span class="title">Sp</span><span class="params">(fvc::div(phiRa)</span></span>, Ua)</span><br></pre></td></tr></table></figure>
<p><code>fvm::laplacian(nuEffa, Ua)</code>对应$\nabla \cdot \left[ \nu_{eff} \nabla U_a \right ]$，<code>fvc::div(Rca)</code>对应$\nabla \cdot \left[ R_{c,a}\right]$。<br><code>phiRa</code>的定义是<br><figure class="highlight stylus"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br></pre></td><td class="code"><pre><span class="line">-fvc::<span class="function"><span class="title">interpolate</span><span class="params">(nuEffa)</span></span>*mesh.<span class="function"><span class="title">magSf</span><span class="params">()</span></span>*fvc::<span class="function"><span class="title">snGrad</span><span class="params">(alpha)</span></span></span><br><span class="line">            /fvc::<span class="function"><span class="title">interpolate</span><span class="params">(alpha + scalar(<span class="number">0.001</span>)</span></span>)</span><br></pre></td></tr></table></figure></p>
<p>相当于$-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}$ <strong>[ <em>注二</em> ]</strong>。<br>于是 <code>fvm::div(phiRa, Ua, &quot;div(phia,Ua)&quot;)</code> 和 <code>fvm::Sp(fvc::div(phiRa), Ua)</code> 便分别对应 $\nabla \cdot\left[ -\nu_{eff} \frac{(\nabla \alpha_a)}{\alpha_a}(\nabla U_a)\right] $ 和 $U_a\left(\nabla\cdot(-\nu_{eff}\frac{\nabla \alpha_a}{\alpha_a}) \right)$。</p>
<figure class="highlight lisp"><table><tr><td class="gutter"><pre><span class="line">1</span><br></pre></td><td class="code"><pre><span class="line"><span class="list">(<span class="keyword">fvc</span>:<span class="keyword">:grad</span><span class="list">(<span class="keyword">alpha</span>)</span>/<span class="list">(<span class="keyword">fvc</span>:<span class="keyword">:average</span><span class="list">(<span class="keyword">alpha</span>)</span> + scalar<span class="list">(<span class="number">0.001</span>)</span>)</span> &amp; Rca)</span></span><br></pre></td></tr></table></figure>
<p>对应$\frac{\nabla(\alpha_a)}{\alpha_a} \cdot \left[ R_{c,a}\right]$，其中<code>&amp;</code>运算符已重载为计算矢量与张量的点乘积 <strong>[ <em>注三</em> ]</strong>。</p>
<figure class="highlight haml"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br></pre></td><td class="code"><pre><span class="line"> =<span class="ruby">=</span><br><span class="line"></span><span class="comment">//  g                          // Buoyancy term transfered to p-equation</span></span><br><span class="line">-<span class="ruby"> <span class="symbol">fvm:</span><span class="symbol">:Sp</span>(beta/rhoa*<span class="constant">K</span>, <span class="constant">Ua</span>)</span><br><span class="line"></span><span class="comment">//+ beta/rhoa*K*Ub             // Explicit drag transfered to p-equation</span></span><br><span class="line">-<span class="ruby"> beta/rhoa*(liftCoeff - <span class="constant">Cvm</span>*rhob*<span class="constant">DDtUb</span>)</span><br><span class="line"></span>        );</span><br></pre></td></tr></table></figure>
<p><code>fvm::Sp(beta/rhoa*K, Ua)</code>对应 $\frac{\alpha_b}{\rho_a} K U_a$，<code>beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)</code> 对应<br>$$<br>\frac{\alpha_b}{\rho_a} \left\{ {C_l (\alpha_b \rho_b + \alpha_a \rho_a) U_r \times (\nabla \times U) -  C_{vm}\rho_b\left[ {\frac{\partial U_b}{\partial t} + U_b \cdot \nabla U_b }\right] } \right\}<br>$$<br>其中变量<code>liftCoeff</code>定义为<br><figure class="highlight css"><table><tr><td class="gutter"><pre><span class="line">1</span><br></pre></td><td class="code"><pre><span class="line"><span class="tag">volVectorField</span> <span class="tag">liftCoeff</span>(<span class="tag">Cl</span>*(<span class="tag">beta</span>*<span class="tag">rhob</span> + <span class="tag">alpha</span>*<span class="tag">rhoa</span>)*(<span class="tag">Ur</span> ^ <span class="rule"><span class="attribute">fvc</span>:<span class="value">:<span class="function">curl</span>(U)))</span></span>;</span><br></pre></td></tr></table></figure></p>
<p><code>DDtUb</code>定义为<br><figure class="highlight aspectj"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br></pre></td><td class="code"><pre><span class="line">DDtUb =</span><br><span class="line">        fvc::ddt(Ub)</span><br><span class="line">      + fvc::div(phib, Ub)</span><br><span class="line">      - fvc::div(phib)*Ub;</span><br></pre></td></tr></table></figure></p>
<p>重力 $g$ 以及曳力的显式项 $\frac{\alpha_b}{\rho_a} K U_b $ 如注释所述，将会在 <code>pEqn</code> 中考虑，压力梯度项 $\frac{\nabla p}{\rho_a}$ 则将在 <code>pEqn</code> 中用来约束两相的连续性。<br>至此动量方程的每一项都与<code>UEqn.H</code>的代码对应起来了。</p>
<h3 id="2-2-_pEqn">2.2. pEqn</h3><p>压力方程的作用是修正两相速度$U_a$ 和 $U_b$以使速度满足连续性方程。将两相的连续性方程加起来，得到总体的连续性方程如下 <strong>[ <em>注四</em> ]</strong>：<br>$$\begin{aligned}<br>&amp; \frac{\partial \alpha_a}{\partial t} + \nabla \cdot (\alpha_a U_a) + \frac{\partial \alpha_b}{\partial t} + \nabla \cdot (\alpha_b U_b) \\<br>= &amp; \frac{\partial (\alpha_a+\alpha_b)}{\partial t} + \nabla \cdot (\alpha_a U_a+\alpha_b U_b) = 0<br>\end{aligned}$$<br>由于$\alpha_a+\alpha_b=1$，于是两相连续性方程等价于<br>$$<br>\nabla \cdot (\alpha_a U_a+\alpha_b U_b) = 0<br>$$</p>
<p>再来看压力方程是如何构建起来的。<br>完整的动量方程离散后，可以写作如下的统一形式：<br>$$<br>a_{p,a}U_{p,a}=H(U_a)-\frac{\nabla p}{\rho_a}+\frac{\alpha_b}{\rho_a} K U_b +g<br>$$</p>
<p>$$<br>a_{p,b}U_{p,b}=H(U_b)-\frac{\nabla p}{\rho_b}+\frac{\alpha_a}{\rho_b} K U_a +g<br>$$<br>其中$H(U_a)$ 和 $H(U_b)$ 包含了动量方程中除 压力梯度项，显式曳力项以及重力项以后所有项的贡献。<br>由此离散方程可以得到 $U_a$ 和 $U_b$ 的表达式如下：<br>$$<br>U_{a}=\frac{1}{a_{p,a}}H(U_a)-\frac{\nabla p}{a_{p,a}\rho_a}+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g<br>$$</p>
<p>$$<br>U_{b}=\frac{1}{a_{p,b}}H(U_b)-\frac{\nabla p}{a_{p,b}\rho_b}+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g<br>$$</p>
<p>如果此 $U_a$ 和 $U_b$ 是方程组的解，那么它们必须满足整体的连续性方程，即<br>$$\begin{aligned}<br>&amp; \nabla \cdot \left[ \alpha_a (\frac{1}{a_{p,a}}H(U_a)-\frac{\nabla p}{a_{p,a}\rho_a}+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g)\right] \\<br>+ &amp; \nabla \cdot \left[ \alpha_b (\frac{1}{a_{p,b}}H(U_b)-\frac{\nabla p}{a_{p,b}\rho_b}+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g ) \right] = 0<br>\end{aligned}$$<br>将压力梯度项移到方程的一边，得到<br>$$\begin{aligned}<br>&amp; \nabla \cdot \left[ \alpha_a (\frac{1}{a_{p,a}}H(U_a)+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g)\right] + \nabla \cdot \left[ \alpha_b (\frac{1}{a_{p,b}}H(U_b)+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g ) \right] \\<br>= &amp;\nabla \cdot \left[ (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b})  \nabla p \right ]<br>\end{aligned}$$<br><strong>这便是压力修正方程的原型</strong>。<br>在<code>pEqn.H</code>中，压力方程其实修正的是界面通量，压力方程迭代收敛以后能保证界面通量的连续性。所以，散度表达式需要根据高斯定理写成界面通量之和的形式：<br>$$\begin{aligned}<br>&amp; \nabla \cdot \left[ \alpha_a (\frac{1}{a_{p,a}}H(U_a)+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g)\right] \\<br>= &amp; (\alpha_a)_f \left[\sum_f(\frac{1}{a_{p,a}})_f H(U_a)\cdot S_f + \sum_f(\frac{\alpha_b K}{ a_{p,a} \rho_a})_f U_b \cdot S_f +\sum_f(\frac{1}{a_{p,a}})_f g \cdot S_f \right ]<br>\end{aligned}$$<br>下标 $_f$ 表示该项将要在代码中用界面上的变量来表示，在OpenFOAM中，即<code>surfaceScalarField</code>，$S_f$ 表示界面的面积矢量，下面的公式里也是一样。<br>$$\begin{aligned}<br>&amp; \nabla \cdot \left[ \alpha_b (\frac{1}{a_{p,b}}H(U_b)+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g)\right] \\<br>= &amp; (\alpha_b)_f \left[\sum_f(\frac{1}{a_{p,b}})_f H(U_b)\cdot S_f + \sum_f(\frac{\alpha_a K}{ a_{p,b} \rho_b})_f U_a \cdot S_f +\sum_f(\frac{1}{a_{p,b}})_f g \cdot S_f \right]<br>\end{aligned}$$</p>
<p>以及<br>$$<br> \nabla \cdot \left[ (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b})  \nabla p \right ] = \sum_f (\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b})_f (\nabla p) \cdot S_f<br>$$<br>下面是<code>pEqn</code>定义了几个跟界面通量有关的变量：<br><figure class="highlight aspectj"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br><span class="line">8</span><br><span class="line">9</span><br><span class="line">10</span><br><span class="line">11</span><br><span class="line">12</span><br><span class="line">13</span><br><span class="line">14</span><br><span class="line">15</span><br><span class="line">16</span><br><span class="line">17</span><br><span class="line">18</span><br><span class="line">19</span><br><span class="line">20</span><br><span class="line">21</span><br><span class="line">22</span><br><span class="line">23</span><br><span class="line">24</span><br><span class="line">25</span><br><span class="line">26</span><br><span class="line">27</span><br><span class="line">28</span><br><span class="line">29</span><br><span class="line">30</span><br><span class="line">31</span><br><span class="line">32</span><br><span class="line">33</span><br><span class="line">34</span><br><span class="line">35</span><br><span class="line">36</span><br><span class="line">37</span><br></pre></td><td class="code"><pre><span class="line">surfaceScalarField alphaf(fvc::interpolate(alpha));</span><br><span class="line">surfaceScalarField betaf(scalar(<span class="number">1</span>) - alphaf);</span><br><span class="line"></span><br><span class="line">volScalarField rUaA(<span class="number">1.0</span>/UaEqn.A());</span><br><span class="line">volScalarField rUbA(<span class="number">1.0</span>/UbEqn.A());</span><br><span class="line"></span><br><span class="line">phia == (fvc::interpolate(Ua) &amp; mesh.Sf());</span><br><span class="line">phib == (fvc::interpolate(Ub) &amp; mesh.Sf());</span><br><span class="line"></span><br><span class="line">rUaAf = fvc::interpolate(rUaA);</span><br><span class="line">surfaceScalarField rUbAf(fvc::interpolate(rUbA));</span><br><span class="line"></span><br><span class="line">Ua = rUaA*UaEqn.H();</span><br><span class="line">Ub = rUbA*UbEqn.H();</span><br><span class="line"></span><br><span class="line">surfaceScalarField phiDraga</span><br><span class="line">(</span><br><span class="line">    fvc::interpolate(beta/rhoa*K*rUaA)*phib + rUaAf*(g &amp; mesh.Sf())</span><br><span class="line">);</span><br><span class="line"></span><br><span class="line"></span><br><span class="line">surfaceScalarField phiDragb</span><br><span class="line">(</span><br><span class="line">    fvc::interpolate(alpha/rhob*K*rUbA)*phia + rUbAf*(g &amp; mesh.Sf())</span><br><span class="line">);</span><br><span class="line"></span><br><span class="line"></span><br><span class="line">phia = (fvc::interpolate(Ua) &amp; mesh.Sf()) + fvc::ddtPhiCorr(rUaA, Ua, phia) + phiDraga;</span><br><span class="line">phib = (fvc::interpolate(Ub) &amp; mesh.Sf()) + fvc::ddtPhiCorr(rUbA, Ub, phib) + phiDragb;</span><br><span class="line"></span><br><span class="line">phi = alphaf*phia + betaf*phib;</span><br><span class="line"></span><br><span class="line">surfaceScalarField Dp</span><br><span class="line">(</span><br><span class="line">    <span class="string">"(rho*(1|A(U)))"</span>,</span><br><span class="line">    alphaf*rUaAf/rhoa + betaf*rUbAf/rhob</span><br><span class="line">);</span><br></pre></td></tr></table></figure></p>
<p><code>phiDraga</code> 和 <code>phiDragb</code> 分别对应 $(\frac{\alpha_b K}{ a_{p,a} \rho_a})_f U_b \cdot S_f +(\frac{1}{a_{p,a}})_f g \cdot S_f$ 和 $(\frac{\alpha_a K}{ a_{p,b} \rho_b})_f U_a \cdot S_f +(\frac{1}{a_{p,b}})_f g \cdot S_f$</p>
<p>由于13-14行的定义，28-29行中的 <code>(fvc::interpolate(Ua) &amp; mesh.Sf())</code> 和 <code>(fvc::interpolate(Ub) &amp; mesh.Sf())</code> 便分别对应的是 $(\frac{1}{a_{p,a}})_f H(U_a)\cdot S_f$ 和 $(\frac{1}{a_{p,b}})_f H(U_b)\cdot S_f$ 。</p>
<p>有了上面的定义，可以看出31行定义的<code>phi=alphaf*phia + betaf*phib</code>便表示了压力方程的左边。</p>
<p>再看33-36行定义的<code>Dp</code>，很显然，表示的是压力方程右边的$(\frac{\alpha_a }{a_{p,a}\rho_a} + \frac{\alpha_b }{a_{p,b}\rho_b})_f$。</p>
<p>有了以上的定义，便可以构建用于修正界面通量的压力方程了：<br><figure class="highlight lisp"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br></pre></td><td class="code"><pre><span class="line">fvScalarMatrix pEqn</span><br><span class="line">  <span class="list">(</span><br><span class="line">     <span class="keyword">fvm</span>:<span class="keyword">:laplacian</span><span class="list">(<span class="keyword">Dp</span>, p)</span> == fvc:<span class="keyword">:div</span><span class="list">(<span class="keyword">phi</span>)</span></span><br><span class="line">  )</span><span class="comment">;</span></span><br></pre></td></tr></table></figure></p>
<p>如上所述，<code>pEqn</code>收敛以后，得到的就是满足连续性的界面通量了，然后再利用求得的界面通量来修正两相的速度，便得到了满足两相连续性的速度：<br><figure class="highlight gherkin"><table><tr><td class="gutter"><pre><span class="line">1</span><br><span class="line">2</span><br><span class="line">3</span><br><span class="line">4</span><br><span class="line">5</span><br><span class="line">6</span><br><span class="line">7</span><br></pre></td><td class="code"><pre><span class="line">Ua += fvc::reconstruct(phiDraga - rUaAf<span class="keyword">*</span>SfGradp/rhoa);</span><br><span class="line">Ua.correctBoundaryConditions();</span><br><span class="line"></span><br><span class="line">Ub += fvc::reconstruct(phiDragb - rUbAf<span class="keyword">*</span>SfGradp/rhob);</span><br><span class="line">Ub.correctBoundaryConditions();</span><br><span class="line"></span><br><span class="line">U = alpha<span class="keyword">*</span>Ua + beta<span class="keyword">*</span>Ub;</span><br></pre></td></tr></table></figure></p>
<p>注意，<code>Ua</code>和<code>Ub</code>为什么是这样来修正呢？回想上面变量定义那个代码段的13-14行，这两行将<code>Ua</code>和<code>Ub</code>分别定义成了$\frac{1}{a_{p,a}} H(U_a)$ 和 $\frac{1}{a_{p,b}} H(U_b)$。<br>回想<code>Ua</code>和<code>Ub</code>的离散方程的统一形式<br>$$<br>U_{a}=\frac{1}{a_{p,a}}H(U_a)-\frac{\nabla p}{a_{p,a}\rho_a}+\frac{\alpha_b}{ a_{p,a} \rho_a} K U_b +\frac{1}{a_{p,a}} g<br>$$</p>
<p>$$<br>U_{b}=\frac{1}{a_{p,b}}H(U_b)-\frac{\nabla p}{a_{p,b}\rho_b}+\frac{\alpha_a}{ a_{p,b} \rho_b} K U_a +\frac{1}{a_{p,b}} g<br>$$</p>
<p>会发现13-14行定义的<code>Ua</code>和<code>Ub</code>都少了几项，所以缺了的这几项的贡献需要在速度修正步骤加回来，而<code>Ua+=</code>后面的<code>fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa)</code>刚好就对应着<code>Ua</code>缺少的那几项。因为经过压力方程修正以后，界面通量是连续的，所以，将缺失的几项对应的界面通量通过<code>reconstruct</code>函数从界面通量重构从对体中心的速度的贡献，便得到了满足连续性的体中心速度了。对<code>Ub</code>也是同样的。</p>
<p>经过以上步骤，便能得到满足整体连续性的两相速度<code>Ua</code> 和 <code>Ub</code>了。</p>
<h2 id="注释">注释</h2><p><strong>注一</strong>：OpenFOAM 里的<code>div</code>函数，字面意义上看起来好像是散度的意思，实际上，<code>div</code>函数执行的是<strong>加和</strong>运算。举例说，对于<code>fvc::div(phia)</code>，<code>phia</code>是<code>surfaceScalarField</code>，其值为<code>(fvc::interpolate(Ua) &amp; mesh.Sf())</code>，即将存储在体中心的<code>Ua</code>插值到每个网格对应的面的面心，然后用面心的速度与该面的面积矢量点乘。从代码中看，<code>fvc::div(phia)</code>对应的是 $\nabla \cdot U_a$，根据高斯定理，也就是 $\sum_f (U_a)_f \cdot S_f$，而<code>phia</code>对应着 $(U_a)_f \cdot S_f$，所以，<code>fvc::div(phia)</code>实际进行的运算是将包围每个网格的面上的通量加起来。更详细的说明见我的<a href="http://xiaopingqiu.github.io/2015/05/17/OpenFOAMcode1/" target="_blank" rel="external">另一篇博文</a>。<br><strong>注二</strong>：注意这里的$\frac{\nabla \alpha_a}{\alpha_a}$ 在代码中的表示方法，详细说明见我的<a href="http://xiaopingqiu.github.io/2015/05/17/OpenFOAMcode1/" target="_blank" rel="external">另一篇博文</a>。<br><strong>注三</strong>：注意这里的$\frac{\nabla \alpha_a}{\alpha_a}$ 在代码中的表示方法，以及与上一个$\frac{\nabla \alpha_a}{\alpha_a}$ 的区别，详细说明见我的<a href="http://xiaopingqiu.github.io/2015/05/17/OpenFOAMcode1/" target="_blank" rel="external">另一篇博文</a>。<br><strong>注四</strong>：这里说的总体的连续性方程指的是<strong>总体体积的守恒</strong>，而不是总体质量的守恒，这二者的差异见Henrik Rusche 的 PHD 论文 P112 的说明。</p>
<h2 id="参考资料">参考资料</h2><ol>
<li>Henrik Rusche， PHD Thesis， Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions， Imperial College of Science, Technology &amp; Medicine, Department of Mechanical Engineering, 2002</li>
<li><a href="https://openfoamwiki.net/index.php/BubbleFoam" target="_blank" rel="external">https://openfoamwiki.net/index.php/BubbleFoam</a></li>
<li><a href="http://www.cfd-online.com/Forums/openfoam-solving/71141-rewriting-twophaseeulerfoam-conservative-form.html" target="_blank" rel="external">http://www.cfd-online.com/Forums/openfoam-solving/71141-rewriting-twophaseeulerfoam-conservative-form.html</a></li>
</ol>

      
    </div>
    <footer class="article-footer">
      <a data-url="http://xiaopingqiu.github.io/2015/05/17/twoPhaseEulerFoam2/" data-id="cj2lay30s000zkkmbrxy8g2n1" class="article-share-link">Share</a>
      
      
  <ul class="article-tag-list"><li class="article-tag-list-item"><a class="article-tag-list-link" href="/tags/Code-Explained/">Code Explained</a></li><li class="article-tag-list-item"><a class="article-tag-list-link" href="/tags/OpenFOAM/">OpenFOAM</a></li></ul>

    </footer>
  </div>
  
    
<nav id="article-nav">
  
    <a href="/2015/05/31/foamTimeAverage/" id="article-nav-newer" class="article-nav-link-wrap">
      <strong class="article-nav-caption">Newer</strong>
      <div class="article-nav-title">
        
          foamTimeAverage
        
      </div>
    </a>
  
  
    <a href="/2015/05/17/twoPhaseEulerFoam1/" id="article-nav-older" class="article-nav-link-wrap">
      <strong class="article-nav-caption">Older</strong>
      <div class="article-nav-title">twoPhaseEulerFoam 全解读之一</div>
    </a>
  
</nav>

  
</article>


<section id="comments">
  <!-- 多说评论框 start -->
  <div class="ds-thread" data-thread-key="post-twoPhaseEulerFoam2" data-title="twoPhaseEulerFoam 全解读之二" data-url="http://xiaopingqiu.github.io/2015/05/17/twoPhaseEulerFoam2/"></div>
  <!-- 多说评论框 end -->
  <!-- 多说公共JS代码 start (一个网页只需插入一次) -->
  <script type="text/javascript">
  var duoshuoQuery = {short_name:"xiaopingqiu"};
	(function() {
		var ds = document.createElement('script');
		ds.type = 'text/javascript';ds.async = true;
		ds.src = (document.location.protocol == 'https:' ? 'https:' : 'http:') + '//static.duoshuo.com/embed.js';
		ds.charset = 'UTF-8';
		(document.getElementsByTagName('head')[0] 
		 || document.getElementsByTagName('body')[0]).appendChild(ds);
	})();
	</script>
  <!-- 多说公共JS代码 end -->
</section>

</section>
        
          <aside id="sidebar">
  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Categories</h3>
    <div class="widget">
      <ul class="category-list"><li class="category-list-item"><a class="category-list-link" href="/categories/C/">C++</a><span class="category-list-count">2</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/DEM/">DEM</a><span class="category-list-count">1</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/OpenFOAM/">OpenFOAM</a><span class="category-list-count">44</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/Paraview/">Paraview</a><span class="category-list-count">5</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/swak4Foam/">swak4Foam</a><span class="category-list-count">1</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/test/">test</a><span class="category-list-count">2</span></li><li class="category-list-item"><a class="category-list-link" href="/categories/vim/">vim</a><span class="category-list-count">1</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Tags</h3>
    <div class="widget">
      <ul class="tag-list"><li class="tag-list-item"><a class="tag-list-link" href="/tags/Boundary-conditions/">Boundary conditions</a><span class="tag-list-count">6</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/C/">C++</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/CentOS/">CentOS</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Code-Explained/">Code Explained</a><span class="tag-list-count">29</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/LES/">LES</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/LIGGGHTS/">LIGGGHTS</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/ODE/">ODE</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/OpenFOAM/">OpenFOAM</a><span class="tag-list-count">20</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Postprocessing/">Postprocessing</a><span class="tag-list-count">9</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Preprocessing/">Preprocessing</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/RTS/">RTS</a><span class="tag-list-count">3</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/TIL/">TIL</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/Windows/">Windows</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/fvOptions/">fvOptions</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/groovyBC/">groovyBC</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/paraview/">paraview</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/test/">test</a><span class="tag-list-count">2</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/thermophysicalModels/">thermophysicalModels</a><span class="tag-list-count">5</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/turbulence-model/">turbulence model</a><span class="tag-list-count">7</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/vim/">vim</a><span class="tag-list-count">1</span></li><li class="tag-list-item"><a class="tag-list-link" href="/tags/wall-functions/">wall functions</a><span class="tag-list-count">4</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Tag Cloud</h3>
    <div class="widget tagcloud">
      <a href="/tags/Boundary-conditions/" style="font-size: 15.56px;">Boundary conditions</a><a href="/tags/C/" style="font-size: 11.11px;">C++</a><a href="/tags/CentOS/" style="font-size: 10px;">CentOS</a><a href="/tags/Code-Explained/" style="font-size: 20px;">Code Explained</a><a href="/tags/LES/" style="font-size: 10px;">LES</a><a href="/tags/LIGGGHTS/" style="font-size: 10px;">LIGGGHTS</a><a href="/tags/ODE/" style="font-size: 10px;">ODE</a><a href="/tags/OpenFOAM/" style="font-size: 18.89px;">OpenFOAM</a><a href="/tags/Postprocessing/" style="font-size: 17.78px;">Postprocessing</a><a href="/tags/Preprocessing/" style="font-size: 11.11px;">Preprocessing</a><a href="/tags/RTS/" style="font-size: 12.22px;">RTS</a><a href="/tags/TIL/" style="font-size: 10px;">TIL</a><a href="/tags/Windows/" style="font-size: 10px;">Windows</a><a href="/tags/fvOptions/" style="font-size: 11.11px;">fvOptions</a><a href="/tags/groovyBC/" style="font-size: 10px;">groovyBC</a><a href="/tags/paraview/" style="font-size: 10px;">paraview</a><a href="/tags/test/" style="font-size: 11.11px;">test</a><a href="/tags/thermophysicalModels/" style="font-size: 14.44px;">thermophysicalModels</a><a href="/tags/turbulence-model/" style="font-size: 16.67px;">turbulence model</a><a href="/tags/vim/" style="font-size: 10px;">vim</a><a href="/tags/wall-functions/" style="font-size: 13.33px;">wall functions</a>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Archives</h3>
    <div class="widget">
      <ul class="archive-list"><li class="archive-list-item"><a class="archive-list-link" href="/archives/2017/05/">五月 2017</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/08/">八月 2016</a><span class="archive-list-count">8</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/06/">六月 2016</a><span class="archive-list-count">5</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/05/">五月 2016</a><span class="archive-list-count">3</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/04/">四月 2016</a><span class="archive-list-count">12</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2016/03/">三月 2016</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/12/">十二月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/11/">十一月 2015</a><span class="archive-list-count">2</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/10/">十月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/09/">九月 2015</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/08/">八月 2015</a><span class="archive-list-count">1</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/06/">六月 2015</a><span class="archive-list-count">2</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/05/">五月 2015</a><span class="archive-list-count">6</span></li><li class="archive-list-item"><a class="archive-list-link" href="/archives/2015/04/">四月 2015</a><span class="archive-list-count">2</span></li></ul>
    </div>
  </div>

  
    
  <div class="widget-wrap">
    <h3 class="widget-title">Recents</h3>
    <div class="widget">
      <ul>
        
          <li>
            <a href="/2017/05/12/duoshuoAnnouncement/">多说评论系统将停止提供服务</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewScritps/">Paraview 脚本一例</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewCamera/">Paraview 中有关 Camera 的操作两例</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewCustomFilter/">Paraview 中创建 Custom Filter</a>
          </li>
        
          <li>
            <a href="/2016/08/27/ParaviewStreamLineOnSlice/">在 Paraview 中画截面上的流线</a>
          </li>
        
      </ul>
    </div>
  </div>

  
</aside>
        
      </div>
      <footer id="footer">
  
  <div class="outer">
    <div id="footer-info" class="inner">
      &copy; 2017 Giskard Q.<br>
      Powered by <a href="http://hexo.io/" target="_blank">Hexo</a>
    </div>
  </div>
</footer>

<script src="//dn-lbstatics.qbox.me/lbservice/busuanzi/2.0/busuanzi.mini.js"/></script>

    </div>
    <nav id="mobile-nav">
  
    <a href="/" class="mobile-nav-link">Home</a>
  
    <a href="/archives" class="mobile-nav-link">Archives</a>
  
    <a href="/atom.xml" class="mobile-nav-link">Rss</a>
  
    <a href="/about" class="mobile-nav-link">About</a>
  
</nav>
    

<script src="//ajax.googleapis.com/ajax/libs/jquery/2.0.3/jquery.min.js"></script>


  <link rel="stylesheet" href="/fancybox/jquery.fancybox.css" type="text/css">
  <script src="/fancybox/jquery.fancybox.pack.js" type="text/javascript"></script>


<script src="/js/script.js" type="text/javascript"></script>

<script>
  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');

  ga('create', 'UA-62501686-1', 'auto');
  ga('send', 'pageview');

</script>

  </div>
<!-- mathjax config similar to math.stackexchange -->

<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    tex2jax: {
      inlineMath: [ ['$','$'], ["\\(","\\)"] ],
      processEscapes: true
    }
  });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {
        skipTags: ['script', 'noscript', 'style', 'textarea', 'pre', 'code']
      }
    });
</script>

<script type="text/x-mathjax-config">
    MathJax.Hub.Queue(function() {
        var all = MathJax.Hub.getAllJax(), i;
        for(i=0; i < all.length; i += 1) {
            all[i].SourceElement().parentNode.className += ' has-jax';
        }
    });
</script>

<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</body>
</html>